GOAT: A simulation code for high-intensity beams 
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A simulation code, GOAT, is developed to simulate single-bunch intensity-dependent effects and their inter- 
play in the proton ring (pRing) of the Electron-Ion Collider in China (EicC) project. GOAT is a scalable and 
portable macroparticle tracking code written in Python and coded by object-oriented programming technology. 
It allows for transverse and longitudinal tracking, including impedance, space charge effect, electron cloud ef- 
fect, and beam-beam interaction. In this paper, physical models and numerical approaches for the four types of 
high-intensity effects, together with the benchmark results obtained through other simulation codes or theories, 
are presented and discussed. In addition, a numerical application of the cross-talk simulation between the beam- 
beam interaction and transverse impedance is shown, and a dipole instability is observed below the respective 
instability threshold. Different mitigation measures implemented in the code are used to suppress the instability. 
The flexibility, completeness, and advancement demonstrate that GOAT is a powerful tool for beam dynamics 
studies in the EicC project or other high-intensity accelerators. 
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1 I. INTRODUCTION 


2 The Electron-Ion Collider in China (EicC) [1] is a pro- 
3 posed highly polarized electron-ion collider based on the 
4 high-intensity heavy-ion accelerator facility (HIAF) [2] up- 
5 grade. It provides a platform for frontier research in nuclear 
e physics with a center of mass energy of 16.7 GeV and lu- 
7 minosity of 2x10°° cm~s~!. The primary parameters are 
s presented in Table 1. The proton ring (pRing), one of the 
9 major accelerators of the EicC project, is both the accelerator 
10 ring and the collider ring of the proton beam. In the pRing, 
11 an Operation mode with multi-bunches and high single-bunch 
12 intensity is essential to achieve the design goal of high lu- 
13 minosity. However, this leads to enormous beam dynamics 
14 Challenges in the pRing. Intensity-dependent effects, such as 
is the beam coupling impedance [3, 4], space charge effect [5], 
16 electron cloud instability [6, 7], beam—beam interaction [8- 
17 10], and their interplay [11-14], are the most severe limita- 
is tions of EicC performance. 

19 To thoroughly investigate the machine performance, beam 
20 instabilities, and associated mitigation methods subjected to 
21 complex high-intensity effects, macroparticle tracking is the 
22 Only way to model and optimize the beam dynamics in pRing. 
23 Various beam dynamics simulation codes have been devel- 
24 oped by CERN according to specific requirements. For ex- 
25 ample, PYHEADTAIL [15] is developed for electron cloud 
26 instability and impedance induced single-bunch collective ef- 
fects, PYECLOUD [16] is exploited for electron cloud es- 
2s tablishment simulation, and the longitudinal collective ef- 
29 fect and beam manipulation are implemented in BLonD [17]. 
30 Beam—beam interactions in colliders can be studied using 
31 BeamBeam3D [18] and Athena [19] developed by LBNL and 
32 IMP, respectively. These codes have been well benchmarked 
33 With different existing codes or beam-based measurements 
34 and have become powerful tools for various beam dynam- 
35 ics studies. However, as the machine performance is con- 
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Table 1. Main parameters of EicC. 


Proton Electron 

Circumference [m] 1341.5928 767.4687 
Kinetic energy [GeV] 19.08 3.5 
Collision frequency [MHz] 100 
Polarization 70% 80% 
Intensity [10' ppb] 1.25 1.70 
Beam current [A] 2.0 2.72 
Bz/B, [m] 0.04/0.02 0.20/0.06 

z 12, im] 10.04/9.58 8.67/7.60 
RMS emittance (H/V) [nmrad] 300/180 60/60 
RMS bunch length [m] 0.04 0.02 
RMS momentum spread [107°] 1.62 0.65 
Transverse tune (H/V) 21.31/22.32 14.08/16.06 
Longitudinal tune 0.0125 0.035 
Laslett tune shift 0.09 - 
Beam-beam parameter 0.004/0.004 0.088/0.048 


Crossing angle [mrad] 50 
Luminosity [em~2s~ 7] 210°" 


3 tinuously pushed to a higher level, different effects can no 
37 longer be treated independently and their complex interplay 
3s Should be considered in any realistic attempt to study the 
39 high-intensity beam dynamics [12, 20]. The combination of 
4 different codes for cross-talk studies is impracticable because 
of differences in coding languages and common interfaces. 
42 Almost no code can simulate multiple high-intensity beam 
43 dynamics in hadron accelerators in a self-consistent man- 
a4 ner. Therefore, as part of the EicC project, a new beam dy- 
4s namics simulation code, GOAT, has been developed in the 
4s past few years and is capable of studying different single- 
47 bunch intensity-dependent beam dynamics, their complex in- 
4s terplay, and possible mitigation techniques simultaneously in 
49 the pRing. 
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5o GOAT is a multiparticle tracking code developed based on 
51 the Python [21] language using object-oriented programming 
52 (OOP) technology. Python significantly improves the effi- 
53 ciency of code development. Some of the core computing 
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modules of the code are written using the Cython [22] pack- 
age because of the slow execution speed of the interpreted 
language. Taking advantage of OOP technology, diversified 
elements and functions can be easily inserted into the code 
without affecting the existing ones. This also facilitates fur- 
ther parallelization to achieve better performance. Thus far, 
based on the most basic beam transportation and manipula- 
tion methods in the transverse and longitudinal planes, im- 
plementations of the impedance induced single-bunch collec- 
tive instability, space charge effect, electron cloud effect, and 
beam-beam interaction are included in the code. 

The remainder of this paper is organized as follows. The 
code architecture and numerical model are given in Section II. 
In Section III to Section VI, the detailed implementations of 
the four different high-intensity effects are explained sepa- 
rately, and the benchmark results against other codes or theo- 
ries are presented and discussed. An application of the code 
for the cross-talk between two different effects is presented 
in Section VII. Section VIII presents the summary and out- 
look. As an important aspect of simulation software, the code 
performance is preliminarily tested and the results are sum- 
marized in Appendix A. 


Il. CODE ARCHITECTURE AND NUMERICAL MODEL 


GOAT consists of a beam module, data management mod- 
ule, infrastructure element module, and physical element 
module. The code architecture is illustrated in Figure 1. All 
accelerator beam particles are defined in the beam module. 
In this module, different particle distributions, such as KV 
and Gaussian distributions, can be generated according to the 
parameters given by the user [23]. Various methods for trans- 
forming the coordinates of the particles and other parameters 
are implemented in the beam module, which can provide suit- 
able parameters for specific studies. In addition, the beam 
module includes functions for the initialization and dynamic 
management of the electrons distributed in the vacuum cham- 
ber used for electron cloud simulation. 

The data management module can read external input files, 
such as beam distribution and optics files according to user 
commands, and convert them into data structures compatible 
with GOAT. The data, such as turn-by-turn particle statistical 
information and instantaneous particle distribution at speci- 
fied time steps during the simulation process, can be stored 
by this module and written to the disk with time stamps. This 
module also includes data post-processing functions, such as 
reading the output files from the simulation, picking up the 
required data, and automatically generating charts. 

Infrastructure element module and physical element mod- 
ule constitute the physical kernel of the GOAT code. The for- 
mer typically contains elements that provide auxiliary func- 
tions for simulations, such as beam slicing elements and PIC 
method-based Poisson solvers. Intensity-independent beam 
transportation and manipulation can also be achieved by com- 
bining infrastructure elements. The latter specifically refers 
to the four types of high-intensity beam dynamics elements 
implemented in GOAT. Each of the four elements is used to 
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Fig. 1. Code architecture of GOAT. 


describe a particular physical process. Based on the OOP 
technology, it is convenient to build new elements in both the 
infrastructure and physical element modules. Instead of using 
a command file to control the simulation procedure [16], the 
user can interactively customize the beamline in a specific 
order by using different combinations of infrastructure and 
physical elements to complete the desired simulation [15, 17]. 

In GOAT, a macroparticle beam of intensity N, energy Eo, 
and particle charge q is described by the 6D set of coordinates 
(£, Px» Y, Py, Z, Pz), where x denotes the horizontal offset 
from the reference orbit, y the vertical offset, p, and p, the 
corresponding transverse normalized momenta, z the longitu- 
dinal offset from the synchronous particle, and p, the relative 
momentum deviation [24]. As for an electron macroparticle 
used in the electron cloud study, a 7D vector (x, Uz, Y, Vy, 
S, Us, Ne) is constructed [6, 16], where x and y are defined 
in the same way as a beam particle, s denotes the distance 
to the transverse reference plane, vz, Vy, and vs the absolute 
velocity component in the horizontal, vertical, and longitu- 
dinal direction, respectively. Further, ne denotes the amount 
of charge carried by the macroparticle. Hereafter, for conve- 
nience, “beam” and “electron” are used to refer to an accel- 
erator beam and electron in the electron cloud study, respec- 
tively. 

Because of the distinctive motion characteristics of the 
beam particles and distributed electrons (the former always 
fulfills ultra-relativistic conditions, whereas the latter is usu- 
ally non-relativistic), different methods are required to inte- 
grate the corresponding equation of motion. GOAT mod- 
els the transverse beam dynamics by tracking macroparticles 
linearly through a transfer matrix between interaction points 
around the ring [24]. Longitudinal beam motion can be mod- 
eled by either linear tracking or nonlinear drift-kick integra- 
tion through RF elements with a given voltage, phase, and 
frequency [25]. However, for an electron, the corresponding 
coordinates are advanced by integrating the Lorentz equation 
that it follows directly. A second-order symplectic Boris algo- 
rithm is implemented in the code [25], and in certain specific 
cases, the integration algorithm can be simplified for better 
performance [7]. 
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In addition, GOAT is equipped with two beam slicing 
methods: one is the uniform-length slicing method, in which 
all slices have the same length, and the other is the uniform- 
charge slicing method, in which each slice contains the same 
number of macroparticles [15]. In each method, slicing is 
performed for the beam rather than the ring, and the upper and 
lower limits of the longitudinal extent are given by the user. 
Both methods have their own advantages and can be used for 
different simulations. Additionally, to meet the requirements 
for solving the Poisson equation, two types of solvers based 
on the particle-in-cell (PIC) method are implemented. One 
uses the finite-difference (FD) time-domain method with per- 
fect electric conducting boundary conditions [6]. The vacuum 
pipe is set as the boundary, which could be either elliptical 
or rectangular. The other method uses the integrated Green 
function (IGF) with open space boundary conditions [9]. The 
boundary is assumed to be at infinity, which is valid when the 
beam size is much smaller than the vacuum chamber. This 
can save computational resources and ensure computational 
speed while maintaining accuracy if the ratio of the pipe size 
to the beam size is large. Furthermore, the linear chromatic- 
ity, the thin nonlinear elements [24], and the bunch-by-bunch 
feedback system [26] are available in GOAT to explore pos- 
sible mitigation measures for high-intensity effects. 

Relying on the abundant infrastructure element module and 
flexible numerical model, the impedance induced collective 
instability, space charge effect, electron cloud effect, beam- 
beam interaction, and their interplay can be simulated in the 
GOAT code using the kick approximation. Together with the 
benchmark results against other frequently-used codes, their 
physical models and numerical approaches are introduced in 
the following sections. 


HI. IMPEDANCE 


In GOAT, the simulation of impedance induced collec- 
tive instability is simple. Only beam transportation and 
impedance elements are required to form the beamline. 
The transverse dipolar impedance, transverse quadrupolar 
impedance, and longitudinal impedance are available in the 
code. Instead of the impedance, the wake function, which is 
the equivalent expression of the impedance in the time do- 
main, is used to calculate the wake kick experienced by the 
beam particles [4, 24]: 


q NSlice 
APu,i BE ` N;[(u;}Wp(zi — Zj) 
j=i+1 
ar uiWo (zi = zi (1) 
ni NStice 
Apz,i B2Eo os Nj; Wo(zi = Zj), 


where u denotes x or y of the ith particle, (u;} the centroid of 
all N; particles in the jth slice, p, the momentum deviation, 
p the Lorentz velocity, and Wp, Wo, and Wo the transverse 
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dipolar wakefield, transverse quadrupolar wakefield, and lon- 
gitudinal wakefield, respectively. 

For the broadband resonator (BBR) impedance model, a 
built-in method is used to calculate the wake function in the 
space domain. For other impedance models, coarse wake 
functions are obtained by reading external input files with two 
columns of data: one for the position and the other for the cor- 
responding wake function. The linear interpolation method 
is adopted for finer wake function calculations. The slicing 
technique can be employed to achieve a better computational 
performance than calculating the wake force between two 
macroparticles [27]. The wake force between slices is calcu- 
lated using the transverse and longitudinal centroids. All the 
macroparticles contained in each slice experienced the same 
wake force. If only the transverse impedance is included, the 
beam is transported linearly in the ring using a transfer ma- 
trix for both the transverse and longitudinal planes. Only an 
RF element is required for the longitudinal impedance study. 
It is noteworthy that the RF element models particle motion 
in the (AE, 0) phase space, whereas the (z, p,) phase space 
is used in the impedance element. The particle coordinates 
are transformed using the beam module. In both cases, the 
wake kick for the beam is integrated at a single interaction 
point [24, 27]. 


A. Transverse mode coupling instability 


In transverse mode coupling instability (TMC), the fre- 
quency shift of each azimuthal satellite increases with 
impedance or beam intensity. When two adjacent modes col- 
lide and merge, the beam becomes unstable, and the oscilla- 
tion of the bunch center starts to grow exponentially [4]. This 
type of instability usually requires a detailed study because it 
is destructive and can cause beam loss. Based on the beam 
parameters in Table 1, the impedance threshold for pRing is 
studied. BBR is used in the simulation to approximate the 
transverse impedance model in pRing, which is a reasonable 
estimation when the pRing’s impedance model is not fully 
built. The following values are considered for other param- 
eters: quality factor Q = 1 and resonant frequency f, = 1 
GHz. The parameters are chosen mainly because the cutoff 
frequnecy of the vacuum chamber with an average radius is 
approximately 1 GHz. 

The proton beam is initialized with transverse and longi- 
tudinal Gaussian distributions. Taking the vertical plane as 
an example, the impedance is scanned from 0 to 8 MQ/m 
at intervals of 0.04 MQ/m. Each set of simulations is per- 
formed for 215 turns. Figure 2(a) shows the spectrum ob- 
tained via fast Fouirer transformation (FFT) using the turn- 
by-turn vertical bunch centroid recorded in the simulation. 
As shown in the figure, with the increase of impedance, the 
frequency of 0 mode and -1 mode moves down and up, re- 
spectively. These two modes are coupled at an impedance 
of 7 MQ/m, after which beam loss occurs rapidly and the 
spectrum is no longer accurate. As a benchmark, the same 
simulation is performed using PyYHEADTAIL, and the spec- 
trum is shown in Figure 2(b). As observed, the behaviors of 
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Fig. 2. Real part of the normalized mode-frequency shifts in the transverse plane given by (a) GOAT and (b) PPHEADTAIL via the macropar- 
ticle tracking with the BBR impedance model. 


200 In addition, in Figure 3, the instability growth rates ob- 
H Parui 261 tained from the two codes are compared by fitting the en- 
GOAT À ee 

f 262 velope of the bunch centroid oscillation. The fitting win- 
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263 dow is in the interval of 20-80% of the total data selected 
264 before beam loss. Again, the TMCI growth rates predicted 
265 by the two codes matches perfectly. All benchmark results 
266 against PPHEADTAIL show that the beam transportation ele- 
267 ment and impedance element implemented in GOAT function 
28 very well. 
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270 In the longitudinal plane, the induced voltage gener- 


271 ated by the longitudinal impedance is superimposed on 
272 the RF voltage, causing potential well distortion. As the 
273 impedance increases further, the microwave instability region 
l aa oN 274 is reached [28]. When the instability threshold is exceeded, 
Fig. 3. Growth rates of the transverse mode coupling instability as a meihe cne spread and banc idene sark aaron dasri 
function of the impedance given by PYHEADTAIL (black dots) and i Ey SP . ; ; 8 ama M 
GOAT (red dots). 276 Although the microwave instability can self-stabilize and will 
277 not lead to beam loss, the beam-beam performance in col- 
278 liding beams is likely to be affected. Similar to the transverse 
279 plane, the BBR is used to estimate the longitudinal impedance 
250 the O and -1 modes are exactly the same as those predicted 20 model in pRing. The quality factor is Q = 1. For the reso- 
251 by GOAT, but the instability threshold of 7 MQ//m is slightly 2s: nant frequency, fy = 2 GHz is chosen such that the rms bunch 
252 higher. The prediction of the instability threshold differs by 2s2 length is longer than the oscillation period of the BBR model, 
253 Only 0.6% between the two codes. Comparing the two spec- 2s3 and makes the bunch particles see mostly the inductive part of 
254 tra, -2 mode and +1 mode are also clearly visible apart from 2s4 the impedance. Thus, the extent of bunch lengthening can be 
255 the two strongest modes. Meanwhile, several lines appeared 285 investigated, as it is an important factor for peak luminosity. 
256 for each azimuthal satellite, which can be interpreted as dif- 2s The simulations are conducted using BLonD and GOAT. 
257 ferent radial modes of the same azimuthal mode. From the 27 The proton beam parameters used in the simulations are listed 


258 perspective of the spectra, the results of the two codes are 2s in Table 1. The beam size, momentum spread, and emittance 
239 are simulated by increasing the longitudinal impedance from 
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0 to 80 kQ with an interval of 5 kQ. Each set of simula- 
tions is performed for 50000 turns. In order to avoid the fil- 
amentation caused by the induced voltage at the initial stage, 
the impedance is elevated adiabatically to maximum within 
the first 8000 turns, which is more than 100 times the syn- 
chrotron period. As shown in Figure 4, the instability thresh- 
old predicted by both codes is 35 kQ. All the results plotted in 
the figures are obtained by averaging the bunch statistics after 
reaching a dynamic balance. In the potential well distortion 
region, the predictions of the longitudinal beam parameters 
from the two codes are identical. Note that the bunch length 
increased and the momentum spread shrank with an increase 
in impedance. This is because of the absence of synchrotron 
radiation in proton beams. Below the instability threshold, 
the evolution of the bunch length and momentum spread must 
ensure that the longitudinal emittance is conserved [29], as 
shown in Figure 4(c). 

However, as the strength of the microwave instability is 
enhanced with higher impedance, divergences gradually ap- 
peared in the simulation results. Taking the momentum 
spread as an example, the maximum relative error between 
different codes is approximately 9%. Similar to results re- 
ported in Ref. [30], bifurcations arise between the different 
simulation codes when the microwave instability become suf- 
ficiently strong. Nevertheless, both BLonD and GOAT pro- 
vides the same instability threshold, and the tendencies of the 
beam parameters obtained through the different codes are not 
significantly different. The correctness of the GOAT code in 
modeling the longitudinal collective effect is verified. In ad- 
dition, other information, such as the bunch distribution and 
detailed evolution of the beam parameters, can also be ex- 
tracted from the simulations. Figure 5(a) shows the longitudi- 
nal bunch shape in the potential well distortion region with an 
impedance value of 25 kQ, and Figure 5(b) shows the evolu- 
tion of the bunch length and momentum spread, where the mi- 
crowave instability has already occurred with an impedance 
value of 45 KQ. 


IV. SPACE CHARGE EFFECT 


The simulation of the space charge effect is straightfor- 
ward. However, the particle motion is continuous under the 
action of the space charge self-field. Applying an integrated 
space charge kick to the beam once per turn leads to incorrect 
results [27]. Therefore, the ring is divided into several seg- 
ments using an external optical input file. These segments can 
be uniform or nonuniform. A space charge element is placed 
at each node. Within each space charge element, the particle 
coordinates are first Lorentz boosted from the co-moving lab- 
oratory frame to the beam frame, and the FD Poisson solver 
is used to compute the electrostatic field slice-by-slice in the 
beam frame. After considering the magnetic field when con- 
verted to the laboratory frame, the space charge kick is ap- 
plied to each particle. The transverse tracking between two 
adjacent space charge elements is linear. When synchrotron 
motion is considered, longitudinal dynamics can be consid- 
ered as either linear or performed by the RF element. 
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Considering the space charge self-field, the envelope equa- 
tions can be written as [31] 


A 2 P KSC 
K 5 x,geo — 0, 
Tr T als)r r3 2s + ry) 
7 3 Kee (2) 
r, + K,(s)r gre = 
Y 4 4 m 2(ro +Ty) 


where r denotes the rms beam widths, K (s) the focusing 
function of lattice, €geo the geometry emittance, €o the per- 
mittivity of vacuum, y the Lorentz factor, and \ the beam 
line density. The coherent quadrupolar tune that describes 
the envelope oscillation can be derived analytically using 
the coupled envelope equations in a smoothly approximated 
ring [31]: 


SC 3 
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where D denotes the coupling parameter, and Qg, are the 
bare tunes in the horizontal and vertical, respectively. 

GOAT code is employed to study the transverse coupling 
phenomenon induced by the space charge effect. As indicated 
in Equation (3), the envelope tunes depend on the bare trans- 
verse tunes of the ring. Thus, seven groups of simulations 
with different initial bare tunes are performed by tracking 
the beam for 512 turns under the action of the space charge 
self-field. The bare vertical tune changes from @,,=21.28 to 
@,=21.34, whereas the horizontal tune is fixed at Q,=21.31. 
The beam is initialized as a coasting beam with a transverse 
KV distribution represented by 1x 10° macroparticles. The 
horizontal and vertical beam sizes are 20 and 10 mm, re- 
spectively. According to numerical convergence studies, 500 
space charge elements are uniformly placed along the ring. 
The smooth approximation is used in the optics file. In the 
simulations, the rms beam size is calculated and recorded at 
each node. Perform FFT on the turn-by-turn beam size at one 
node, and two envelope coherent modes are obtained. In Fig- 
ure 6, the simulated coherent tunes are compared with those 
predicted by Equation (3). The variation in the vertical beam 
size owing to the vertical beta function scaling with the ver- 
tical tune is considered for both the simulation and analytical 
calculation. As shown in the figure, GOAT matches the the- 
ory well. 

In addition, the incoherent tune spread of the beam parti- 
cles with different amplitudes is a direct dynamic result of 
the space charge self-field. It can be used as a crosscheck for 
space charged elements. Under the smooth approximation, 
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Fig. 4. (a) Bunch length, (b) momentum spread, and (c) emittance as a function of the impedance obtained by GOAT (black dots) and BLonD 
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Fig. 5. (a) Longitudinal bunch shape normalized to unit with the longitudinal impedance value of 25 kQ in the PWD region. (b) Evolution of 
the bunch length and momentum spread with the impedance value of 45 kQ in the microwave instability region. 


all particles in the KV distribution have exactly the same in- 
coherent tune shift when subjected to the space charge self- 
field, which can be expressed by the following formula [31]: 


KSC R? 
Rey(Re T Ry)Qz,y ? 


where Rg, y are the rms beam sizes of a KV distribution in 
horizontal and vertical plane, respectively. Here, a KV beam 
with a radius of 10 mm is used, and the incoherent tune shift 
calculated using Equation (4) is —0.04. Correspondingly, 
the maximum incoherent tune spread of a Gaussian beam can 
also be obtained by analytical calculations [31]: 
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where £ (s) is the beta function, the cx, the rms beam size 
of a Gaussian distribution in the horizontal and vertical plane, 
respectively. According to the theory in accordance with the 
rms equivalence principle, a Gaussian beam with an rms size 
of 5 mm is considered, and the maximum incoherent tune 
spread is —0.08, which is twice the incoherent tune shift in 
the KV beam [5, 32]. For comparison, the numerical track- 
ing is performed using the beam parameters employed in the 
analytical evaluation. In both cases, the beams are initial- 
ized with 5x10° macroparticles, and 100 space charge ele- 
ments are uniformly distributed along the ring. The beams 
are tracked linearly in the transverse plane and the longitu- 
dinal motion is frozen. As shown in Figure 7, the black and 
red dots represent the incoherent tune spreads of all particles 
in the Gaussian and KV beams, respectively. The incoherent 
tune spread of the Gaussian beam particles reaches the the- 
oretically calculated value indicated by the blue pentagram. 
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Fig. 6. Envelop mode tunes as a function of the bare vertical tune of 
the machine for the coasting pRing beam obtained by theory (lines) 
and macroparticle tracking (dots). 
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Fig. 7. Incoherent tune spread obtained by the macroparticle track- 
ing with transversely a Gaussian distribution (black dots) and a KV 
distribution (red dots). The incoherent tune shift of the KV distribu- 
tion is denoted by the green pentagram, and the maximum incoher- 
ent tune spread of the Gaussian distribution is marked by the blue 
pentagram. The bare tune is represented by the purple pentagram. 


Although the tune spread of the KV beam given by the sim- 
ulation is not strictly a point marked by the green pentagram, 
as predicted by theory, the dispersion is very small. This is 
consistent with the results reported in Ref. [5]. This may be 
caused by fluctuations in the particle distribution or numeri- 
cal errors in the field solver. Overall, the space charge element 
implemented in GOAT is correct. 
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Table 2. Numerical parameters used for electron cloud simulations. 


Bunch intensity, Np [1017] 1.25 
Beam size, Ox, [mm] 3.0/2.3 
Chamber shape elliptical 
Horizontal chamber aperture, a [mm] 45 
Vertical chamber aperture, b [mm] 35 
Dipole field [T] 1.0 
Quadrupole gradient [T/m] 20.0 
Electron reflectivity Ro 0.7 
E° for elastic scatters [eV] 150 
s for true secondaries 1.35 
Energy of dmaz, Emax [eV] 332 
Maximum SEY, dmaz 2.0 


V. ELECTRON CLOUD EFFECT 


In modeling electron cloud effects, the cloud build-up pro- 
cess and beam-cloud interaction are treated separately, or 
the so-called “weak-strong” model [6]. This is because the 
dynamic balance of the electron cloud density in the vac- 
uum chamber can be reached in a single or very few turns, 
whereas the beam-cloud instability study requires tracking 
the beam for many turns, at least longer than the instabil- 
ity growth time. In GOAT, two methods are established for 
beam-cloud interactions: the linearized method [31] and the 
self-consistent tracking method [7]. The linearized method is 
based on cloud dipolar and quadrupolar forces obtained from 
dedicated simulations. In the self-consistent tracking method, 
both the beam particles and electrons are characterized by 
macroparticles. The electromagnetic interaction between the 
beam and cloud is modeled by a set of thin interactions along 
the ring, and the forces acting on each other are solved numer- 
ically using the FD Poisson solver implemented in the code. 
Currently, the linearized method is preferred for studying the 
beam-cloud instability because the computation overhead for 
the tracking method is expensive. 


A. Build-up simulation 


In the build-up simulation, the very low density cold elec- 
trons are uniformly generated in the chamber as the “seed” 
electron at the initialization stage. When the bunch passes, 
primary electrons are created due to the residual-gas ioniza- 
tion or beam loss. These electrons, along with the surviving 
electrons, are accelerated to some energy under the action of 
the beam field. Secondary electrons are produced as ener- 
getic electrons hit the chamber wall. Only elastic scattering 
and real secondary emission are included in the secondary 
electron emission model. The energy carried by the electron 
determines the secondary electron yield (SEY) of an elasti- 
cally scattered electron [6]: 
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where Ro corresponds to the elastic reflection probability for 
an electron in the limit of zero primary electron energy, E65 
the fitting parameter of experimentally measured energy spec- 
trum of elastically scattered electrons (chosen as a specific 
value according to the vacuum chamber’s material), and Æ 
the energy of the incident electrons (calculated via the veloc- 
ity of each electron in the simulation, in unit of eV), or a true 
secondary electron is produced with SEY [6]: 
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where Emax is the energy of the electrons at which the SEY 
is Maximum, maz the maximum value of the SEY at nor- 
mal incidence, and s the fitting parameter of measured SEY 
curve (s=1.35 is used widely since it provides the most rea- 
sonable fit to the experimental data.). The newly generated 
electrons follow a logarithmic normal distribution in energy 
and a cosine distribution in direction. This process is repeated 
successively according to the user-defined beam filling pat- 
tern. It is worth noting that, due to the multipacting effect, the 
number of electron macroparticles and the charge carried by a 
single electron change constantly during the simulation [16]. 
For considerations of computing performance and memory, 
as in PYECLOUD, both the number of macroparticles and the 
amount of charge of single macroparticles are dynamically 
managed. All alive electrons are mixed and redistributed in 
the 6D phase space when the above two quantities exceed 
the pre-defined threshold. After the regeneration process, the 
charge of each electron is at the same level. 

Electron cloud build-up simulations are performed in the 
drift, dipole, and quadrupole regions using PyYECLOUD and 
GOAT software. The numerical parameters in Table 2 are 
used. Figure 8 shows the evolution of the electron line den- 
sity within the vacuum chamber simulated by the two codes. 
The line density in the quadrupole region shown in the figure 
is halved for comparison. Both codes yield the same results. 
The time scale shown here is only 1/8 of the revolution pe- 
riod. Nevertheless, owing to the extremely short bunch spac- 
ing, the electron cloud reaches equilibrium after dozens of 
bunches pass through. As the cloud becomes saturated, the 
electron densities in the field-free and quadrupole regions are 
relatively close, at approximately 4x10'°/m. Although the 
horizontal motion is frozen in the presence of a dipole mag- 
netic field, the electron density also reaches 2.5x 101° /m and 
tends to increase slowly. In addition, the oscillation ampli- 
tude of the electron density is smaller in the presence of an 
external magnetic field than in the field-free region because 
most electrons are trapped by the magnetic field lines. 

Figure 9 shows the electron spatial distribution exported 
from GOAT in the dipole region just before the bunch ar- 
rival. Similar to other studies, vertical stripes are formed 
in the horizontal direction in the presence of a dipolar mag- 
netic field [6, 7]. In Figure 10, comparisons of the transverse 
phase space and corresponding histograms of the two codes 
are shown. In the horizontal direction (Figure 10(a)), the elec- 
tron distribution is figure-8 shaped. The central density is rel- 
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Fig. 8. Comparisons of the evolution of the electron line density in 
the drift, dipole, and quadrupole regions. The results are obtained 
by PyECLOUD (solid lines) and GOAT (dashed lines). The electron 
line density in the quadrupole region is halved. 
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Fig. 9. Snapshot of the electron distribution in the dipole region 
before the bunch arrival at the cloud sturation stage. 


atively low. Almost all the electrons gather in the range of two 
stripes. The velocity distribution of the cloud is Gaussian- 
like, and the electrons are nearly static. In the vertical di- 
rection (Figure 10(b)), however, the electrons are distributed 
over the chamber since the motion is unconstrained. At the 
same time, it can be noted that the electrons are highly con- 
centrated near the chamber wall, which is a significant fea- 
ture of a saturated cloud. The electrons in the vicinity of the 
chamber wall form a potential barrier to prevent the energetic 
electrons from hitting the pipe and producing secondary elec- 
trons. The vertical velocity component is sharper than the 
horizontal one, and the value is approximately an order of 
magnitude higher. The comparison shown in the histograms 
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Fig. 10. The horizontal (a) and vertical (b) phase space distributions and statistical histograms of the electron cloud in the presence of a dipole 


magnetic field. The histograms are normalized to unit. 


further confirms the consistency between PyYECLOUD and 
GOAT. 


B. Beam-cloud interaction 


Two quantities are required due to the electron pinch as the 
bunch passes through in modeling the beam-cloud interac- 
tion using the linearized method: the generalized 2D dipolar 
wakefield [33-35] and the longitudinally varied betatron tune 
shift [32, 36]. The dipolar term is first dicussed. When an on- 
axis bunch passes through a cloud, the electrons are attracted 
toward the axis and the centroid of the cloud does not move. 
However, the situation changes if a small transverse offset is 
introduced to the bunch. When the displaced bunch passes 
through a cloud, the cloud is redistributed and begins to os- 
cillate. Then, the subsequent portions of the bunches are de- 
flected. Therefore, to obtain the dipolar wakefield, the bunch 
is sliced longitudinally, and a transverse offset is added to the 
driving slice. The dipolar wakefield generated by the driving 
slice can be computed for the subsequent testing slices via 
[35] 


Ej xy L 


ix,yt\i 


(8) 


where L denotes the length of the cloud, Ey y the electric field 
averaged over the particle distribution, A the transverse offset 
attached to the driving slice containing N particles, and the 
subscript 7 and 7 are respectively the driving and the testing 
slices. Instead of depending on the relative distance between 
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the ¿th driving slice and the jth testing slice only, the dipo- 
lar wakefield produced by the electron cloud depends on two 
longitudinal positions simultaneously as a result of the elec- 
tron pinch. For example, in the field-free region, by chang- 
ing the position of the driving slice along the slice sequence, 
the generalized dipolar wakefield can be obtained, as shown 
in Figure 11. z > 0 corresponds to the bunch head. The elec- 
tron distribution used in the simulation is extracted from the 
dedicated buildup simulation described in the previous sec- 
tion. The transverse offset is 10% of the vertical beam size. 
Longitudinally, a Gaussian distribution with nominal bunch 
length is used. The bunch is cut into 50 slices within the 
range of (—40,, 40), where oz is the rms bunch length. It 
is worth noting that for the purpose of reducing the numerical 
noise in the simulation, the electron distribution is created by 
merging ten realistic distributions at the same time step just 
before the bunch arrival. Moreover, a set of reference “driving 
slice” groups without attached initial transverse offset is es- 
tablished in each simulation, which has the same distribution 
as the beam and is used as the base of the numerical noise. 


As shown in Figure 11, the consequences of the elec- 
tron pinch are significant. The generalized dipolar wakefield 
does not satisfy the translational invariance. This is a ba- 
sic property of the conventional electromagnetic impedance. 
The generalized wakefield peaks at the bunch head at a value 
of 8x10!° VC~!m~!. In contrast to the results reported in 
[35] with an initially static electron distribution, no signifi- 
cant increase in the wakefield amplitude caused by the elec- 
tron pinch is observed when a realistic electron distribution 
is used. This illustrates that the enlargement of the wakefield 
amplitude is closely related to the electron velocity distribu- 
tion. 
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Fig. 11. The generalized 2D dipolar wakefield generated by the elec- 
tron cloud in the field-free region. z > 0 means the head of the 
bunch. The realistic distribution of electrons from dedicated build- 
up simulation are used. 


To obtain the betatron tune shift variation along the bunch 
length, the self-consistent macroparticle tracking method is 
used. A symplectic expression has been derived for the beam- 
cloud interaction [37] 
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where ® is the electrostatic potential of the cloud. Instead of 
the typical approaches of including only the transverse part of 
this map, the complete map is implemeted in GOAT. Based 
on the numerical convergence studies, 20 thin electron cloud 
elements are uniformly placed around the ring. The beam 
is represented by 5x 10° macroparticles with 50 longitudinal 
bunch slices. Fresh cloud distributions with approximately 
2x10° macroparticles at each node are generated and saved 
at the initialization stage. Again, a realistic electron distribu- 
tion is used. Between electron cloud elements, the beam is 
linearly transported in the transverse plane. And the longitu- 
dinal motion is frozen to obtain the instantaneous tune spread. 

In Figure 12, the blue dots represent the incoherent tune 
spread caused by the electron cloud force. There are several 
steps in the tune spread distribution, particularly at the bunch 
head with z > 0. This is mainly because the discrete equa- 
tion of motion is used to integrate the electron motion, and 
the change in the spatial morphology of the cloud is discon- 
tinuous. The tune shift obtained by averaging over the tune 


603 


604 


605 


606 


607 


608 


609 


610 


61 


=e 


612 


613 


614 


615 


616 


617 


618 


619 


620 


62 


mre 


622 
623 


624 


625 


626 


6 


N 


7 


628 


10 


spread of particles within the slices at different longitudinal 
positions is plotted as the red line. The electron pinch effect 
results in a change in the focusing forces during the bunch 
passage, and therefore, in the betatron tune modulation with 
longitudinal coordinates. One can observe that the tune shift 
increases by a factor of eight from the head to the center of 
the bunch and then slowly decays toward the tail. In addition, 
a line parallel to the horizontal axis at AQ,,=0.18 is clearly 
visible. This is because some particles at the bunch center 
experiencing very strong cloud-gradient forces and crossing 
the half-integer resonance Q,,+0.18=22.5, as indicated by the 
dashed green line. 


Taking the above results as a wakefield file for the ex- 
ternal input, the impedance element can be used to simu- 
late the beam-cloud instability. The cloud-generated dipolar 
and quadrupolar forces can be modeled by transverse driving 
and detuning impedance implemented in the impedance ele- 
ment. The impedance element is well benchmarked in Sec- 
tion III. The nonlinear characteristics of the beam-cloud in- 
teraction are omitted in such a linearized method. However, 
this method is correct and has advantages obviously in fast 
and enormous parameter scans [32]. 


VI. BEAM-BEAM INTERACTION 


For colliding beams, the GOAT code describes the electro- 
magnetic interaction of two counter-rotating beams through 
the 6D symplectic Synchro-Betatron Mapping method [38]: 
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Fig. 12. Incoherent tune spread (blue dots) of the bunch under the 
action of the electron cloud force in the field-free region. The varied 
transversal tune shift in y direction (red line) along the longitudinal 
position can be obtained by averaging the incoherent tune spread of 
particles in the slices at each longitudinal position. The half-integer 
resonance line Q,+0.18=22.5 is shown in green dashed line. 
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where the S denotes the drift distance, BỌ P, EOP, and EGP 
the electric field generated by the opposite bunch slice at the 
collision point. The linear transfer matrix is used for the 
633 motion of beam particles in both transverse and longitudinal 
634 planes between different interaction points. 

The calculation sequence of the beam-beam interaction is 
as follows. The two bunches are first sliced longitudinally 
using the uniform charge slicing method implemented in the 
code, and then the bunches collide slice-by-slice. At each 
collision step, the particles contained in the slice are drifted 
from the interaction point (IP) to the collision point (CP) lo- 
cated at S = (z1 + 22) /2, where zı and z are the longi- 
tudinal centroids of the two slices. The IGF Poisson solver 
is used to compute the electric fields generated by the slices 
according to the current distributions. Subsequently, taking 
into account the magnetic field and time-of-flight effects, the 
kick is applied to the particles in the opposite slice. Finally, 
the particles are transferred back to the IP via the inverse drift 
operator. In general, the linear interpolation method [39] is 
applied to calculate the kick for each particle in terms of com- 
putational convergence and speed. 

In addition, the requirements for the fast separation of the 
two colliding beams, the overall detector component and in- 
teraction region (IR) magnet arrangements strongly depend 
on a large crossing angle shown in Figure 13. To include the 
e55 Crossing angle in the beam-beam simulation procedure given 
656 by Equation (10), which is derived without the crossing angle, 
es7 the Lorentz boost [40]: 
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es9 can be used to transform the beam particles from laboratory 
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eo frame to head-on frame. A variable marked with and without 
1 a superscript tilde indicates the quantity in the head-on frame 
and laboratory frame, respectively. The variable 0, denotes 
the half crossing angle, and the h hamiltonian, 
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e5 After collision, the beam particles need to be transformed 
ess back to the laboratory frame via inverse Lorentz boost, 


De = Pz COSA, + h tan 0e, 
Py = Py COs be, 


Dz = Dz + Pz tan be — h tan? be, 
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ess Currently, only the ‘strong-strong’ model is available in 
e9 GOAT. 

670 Predicting beam parameters such as luminosity, beam size, 
671 and centroid oscillation, studying the dependencies between 
672 different parameters, and investigating the role of beam-beam 
673 dynamics on the cross-talk of multiple beam dynamics are 
674 the main tasks of the beam-beam simulation. In order to 
67s Check the validity of the beam-beam model implemented in 
67s GOAT, a set of beam-beam simulations is carried out with 
677 and without considering the crossing angle utilizing Athena 
and GOAT. The beam parameters in Table 1 are used. The 
luminosity and beam size are shown in Figure 14(a) and Fig- 
ure 15. When the two beams collide in the head-on frame, 
the transverse beam sizes of the proton beam are stable at 
the design value due to its small beam-beam parameter. For 
the electron beam, the beam sizes in both the horizontal and 
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Fig. 13. Beams colliding at an angle in the original laboratory frame 
(top) and in the boosted head-on frame (bottom). 


vertical directions shrink to slightly smaller values than the 
design value, which is also why the luminosity is higher than 
the calculated value. For the collision with a crossing angle 
of 50 mrad, although the luminosity is reduced by more than 
six times compared to the head-on collision, the beams are 
stable, and no significant beam sizes blow up are observed. 
The transverse beam sizes of protons and electrons are sta- 
ble around the design values after reaching equilibrium. This 
indicates that the luminosity degradation is caused by geo- 
metric loss. Again, the simulation results given by Athena 
and GOAT agree very well with each other. 


In the EicC conceptual design, the luminosity reduction 
caused by the crossing angle is compensated by the crab cav- 
ity. The transfer map of the crab cavity in the x-z plane for 
each particle is given by [41], 
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Fig. 14. (a) The luminosity evolution in case of collision with 
and without considering the crossing angle obtained by Athena and 
GOAT. (b) The impact of the crab cavity frequency on the luminos- 
ity obtained by GOAT. 


old 
po — old 4 q sin We ate db 
x x BE P 
ld qVw old wer a 
new _ „o o 
Pz = pP; + B2cEo Os ( c +6), 


where V is the voltage of the crab cavity, w the angular fre- 
quency of the cavity, and @ the phase of the cavity. The cavity 
voltage is chosen as: 


cE tan (2s) 


V= 
quy/B* Bco sin Aw’ 


where ĝ* is the beta function at the IP, Bcc the beta func- 
tion at cavity location, and Aw the phase advance between 
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Fig. 15. The corresponding horizontal (a) and vertical (b) beam sizes 
of colliding beams for the cases of Fig. 13(a). 


Fig. 16. The x-z distributions of the bunch at IP with different crab 
cavity frequencies. The axis are normalized to the nominal horizon- 
tal beam size and bunch length. The crab frequency is 100 MHz 
(top left), 200 MHz (top right), 300 MHz (bottom left), and 400 
MHz (bottom right). 
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the IP and the cavity. For the proton beam, the higher en- 
ergy and smaller beta function at the IP lead to a much higher 
cavity voltage requirement than that for the electron beam. 
Generally, higher frequencies are preferred for reducing the 
cavity voltage. Figure 14(b) shows the luminosity evolution 
obtained by scanning the crab cavity frequencies at a fixed 
half-crossing angle of 25 mrad. The luminosity degradation 
is low since the cavity frequency is below 200 MHz. As the 
frequency continuously increases, the luminosity degradation 
becomes higher. This phenomenon can be explained by the 
x-z distribution at the IP for different crab cavity frequencies, 
as shown in Figure 16. Because of the incomplete deflec- 
tion caused by the crab with a higher frequency (lower wave- 
length), the head and tail of the bunch are farther off from the 
axis at the IP, resulting in an increase in the overlap area of 
the two bunches and a decrease in luminosity. In addition, the 
simulated luminosity over 20000 turns is higher than the de- 
sign value, even for a cavity frequency of 400 MHz. However, 
the synchrotron-betatron (SB) resonance may still be excited 
owing to the longitudinally position-dependent beam-beam 
kick, causing a reduction in the luminosity lifetime [42]. 


VII. APPLICATIITON OF CROSS-TALK 


The physical elements implemented in GOAT are well 
benchmarked in the previous subsections. Thus, its validity 
is guaranteed. However, the beam in a real machine cannot 
be affected by a single effect, and there is cross-talk between 
different effects. Among the many high-intensity effects in a 
collider, the beam-beam interaction is one of the most impor- 
tant beam dynamics processes, which has a direct impact on 
the luminosity and integrated luminosity of the machine. Sev- 
eral instabilities have been observed in previous studies due 
to the cross-talk between the beam-beam interaction and other 
intensity-dependent effects [11—13, 43, 44]. It is necessary to 
explore the mechanisms of these instabilities and correspond- 
ing mitigation measures at the machine design stage. 

In this subsection, a numerical example is presented for 
the cross-talk between the beam-beam interaction and the 
pRing’s vertical impedance in EicC. To emphasize the impor- 
tance of cross-talk simulation, three sets of simulations are 
performed. The first purely includes the effect of the pRing’s 
vertical impedance, the second considers only the beam-beam 
interaction, and the third considers the self-consistent treat- 
ment of the two effects. The first two sets of simulations 
are performed in previous subsections, and the results of the 
evolution of beam parameters are presented here. In the 
self-consistent simulation, the transverse impedance element, 
beam-beam interaction element, and transverse and longitu- 
dinal linear transportation elements are combined to form a 
ring. A vertical impedance value of 6 MQ)/m is chosen, 
which is below the instability threshold. And the beam-beam 
interaction is assumed to take place in the head-on frame. 
Figure 17(a) and Figure 17(b) show the evolution of the beam 
vertical centroid and emittance for the three cases. It can be 
observed that below the instability thresholds of TMCI and 
beam-beam interaction considered alone, a dipole instability 
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Fig. 17. Evolution of vertical bunch centroid (a) and emittance (b). 
Three cases are considered, first is purely beam-beam interaction, 
second is purely transverse impedance, and third is cross-talk of 
beam-beam interaction and transverse impedance. 


arises in the simulation. The vertical centroid grows expo- 
nentially, and the emittance blows up. The simulation results 
suggest that there is no coupling between the coherent beam- 
beam mode and the longitudinal sideband because the bare 
tunes of the two beams are different. This is very different 
from the mode coupling instability for symmetric collisions 
reported in Ref. [12]. It is also worth noting that such a dipole 
instability does not appear in the horizontal plane. This im- 
plies that the obseved head-tail type instability is closely re- 
lated to the hourglass effect, because it is the only major dif- 
ference between the two transverse planes. Further studies 
are still required to identify and understand the underlying 
mechanisms of this coherent instability. 


Obviously, this numerical application indicates the ne- 
cessity of simulating multiple physical processes in a self- 
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Fig. 18. The impact of chromaticity on the growth rate of coherent 
instability. 


consistent manner. This also demonstrates the flexibility, 
completeness, and advancement of GOAT. 

In addition, the linear chromaticity and the ideal bunch-by- 
bunch feedback system implemented in GOAT are utilized to 
suppress this coherent instability. In Figure 18, the instability 
growth rates for different chromaticity values are presented 
in the range from -10 to 10. The growth rates increase and 
then remain almost constant when the chromaticity values are 
negative. For positive chromaticity values, non-monotonic 
growth rate behavior is observed. When the charomaticity 
value is greater than eight units, the instability is fully sup- 
pressed. Then, the ideal bunch-by-bunch feedback system 
is employed to damp the instability. The evolution of the 
vertical centroid and emittance for different damping rates 
are presented in Figure 19(a) and Figure 19(b), respectively. 
Compared to chromaticity, the feedback system is more ef- 
fective in suppressing this instability. It can be seen that the 
dipole motion is suppressed, and the emittance is conserved, 
even though the damping rate is very small. This can be ex- 
plained by Figure 20, which illustrates the intra-bunch motion 
of the proton beam extracted from the simulation in succes- 
sive turns. The most unstable mode is the 0-mode and the 
bunch-by-bunch feedback can eliminate the oscillations of the 
bunch as a whole. 


VII. SUMMARY AND OUTLOOK 


In this paper, a simulation code, GOAT, is developed for 
single-bunch high-intensity beam dynamics. It can be used 
to simulate all intensity-dependent effects in the pRing of the 
EicC project. The code architecture and numerical model are 
introduced. Four simulation examples, including impedance 
induced collective instability, space charge effect, electron 
cloud effect, and beam-beam interaction, are conducted based 
on abundant elements and the flexible numerical models pro- 
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Fig. 19. The impact of ideal bunch-by-bunch feedback system on 
the coherent instability. The evolution of vertical bunch centroid (a) 
and emittance (b) with different damping rates. 


vided by GOAT. The results are well benchmarked with other 
codes and theories. In addition to separate simulations, an 
application is presented for cross-talk between the beam cou- 
812 pling impedance and the beam-beam interaction. The com- 
sis prehensiveness and correctness of GOAT are verified. The 
814 Python program coded by the OPP technology ensures the 
sis scalability of GOAT as well as the independence of the mod- 
s16 ules. Different effects can easily be integrated into the code. 
si7 GOAT can also be used to simulate the intensity-dependent 
sis beam dynamics in other accelerators and colliders. In the fu- 
si9 ture, it is planned to improve the performance of the GOAT 
s2 code, including algorithm optimization and hardware support 
s21 (such as parallelization based on GPU (Graphic Processing 
s22 Unit)). With the help of parallelization techniques, GOAT 
ses Can be upgraded to include the multi-bunch beam dynamics. 
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Fig. 20. Intra-bunch motion for successive 100 turns of proton beam 
extracted from the cross-talk simulation. 
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Appendix A: Estimation of the code performance 


The code performance is a crucial aspect of simulation 
software. Due to the abundant physical modules implemented 
in GOAT, a substantial number of numerical parameters con- 
tribute to the computational speed. Therefore, three typical 
quantities, namely the number of macroparticles, slices, and 
mesh grids in the field solver, are selected as test quanti- 
ties. Table 3 summarizes the performance of all modules men- 
tioned in the paper. All tests are run on Intel(R) Core(TM) 
19-9900K CPU @ 3.60GHz (RAM 16.0 GB) adopting single- 
core and single-thread. Some notes are attached to the test 
results. (1) The slicing is performed on the ring in the simula- 
tions of space charge effect and electron cloud buildup, while 
on the beam for the other three simulations. (2) The saturated 
electron cloud is usually established after the passage of mul- 
tiple bunches within one revolution period, thus the time spent 
for each bunch passage is used to describe the elapsed time of 
the electron cloud simulation. Also, it is more intuitive and 
convenient to use the time step corresponding to each slice 
to describe the number of slices used in the electron cloud 
effect simulation. (3) As mentioned, in electron cloud simu- 
lations, the number of electron macroparticles and the amount 
of charge carried by each macroparticle are constantly chang- 
ing due to secondary electron production, and there is no fixed 
number of macroparticles. Two typical preset values given 
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Table 3. The computation performance of each physical module under different numerical parameters. 


(a) macroparticle, Nmp (b) Slice, Ns (c) Grid in field solver, Ng 
Number Time** (s) Number Time** (s) Number Time** (s) 
1x10° 0.0156 10 0.107 2 = 
TMCI 5x10°* 0.109 25* 0.109 = 2 
1x10° 0.237 50 0.111 > 2 
1x10° 0.0313 25 0.188 - - 
LMWI 5x10°* 0.189 50* 0.189 = z 
1x108 0.392 100 0.195 : - 
1x107 0.509 5x5* 1.514 128x128 1.328 
BB 5x105* 1.514 7x7 2.157 256x256* 1.514 
1x10° 3.513 9x9 2.844 512x512 2.974 
1x10° 2.2 50 6.1 128x128 11.4 
SC 5x105* 11.9 100* 11.9 256x256* 11.9 
1x10° 23.7 300 35.6 512x512 15.6 
EC 2.5-5x 107 52 dt=11.25 ps* 68 256x256 37 
5-10x 10** 68 dt=22.5 ps 45 512x512* 68 


* The variable is fixed when the other types of numerical parameters are varying. 
** The elapsed time refers to [s/turn] for TMCI, LMWI, BB, and SC simulations, and [s/bunch] for EC simulation. 


by the user during the initialization stage are used. The first 
number indicates that: when the number of macroparticles ex- 
ceeds this value, the coordinates and charges of all particles 
are mixed and regenerated, and the meaning of the second 
number is: the total number of macroparticles after regenera- 
tion. 


For impedance induced collective effects, the main limi- 
tation comes from the number of macroparticles; for other 
effects, the increase in all three test quantities significantly 
increases the computation time. 
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